/***************************************************************
                      cb_attic.c
was called zutil.c

Zrule_menu it the attice menu
Zgram writes a gram file
CTG_resources - gettime settime
Distribution code written by Jhat, access from main analysis attic
***************************************************************/
#include "clam.h"
#include <fcntl.h>
#include <pwd.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <float.h>

#define SEEK_SET 0
#define BUFSZ 120
extern void getdatehead();
extern char projectname[80];
extern BOOL update;
extern void histZap(), Zsel_unbury(), Zall_unbury(), Zsel_bury(), Zsel_allbury();
extern void C1zBur(), Zexpert_zap(), Select_and_Colour();
extern int loadCzfast();
extern int Zget_tol();
extern void ClamCtgChr();
extern void freelistmem();
extern FILE *graphQueryOpen();

extern void graphArray(Graph *g, int *xarray, int *yarray, int arraysize,
		       char *xtitle, char *ytitle, char *fn, char *title,
		       void (*pick)(int));
extern void setDisplayType(int);
extern void displayProj();
extern void updateproj();
void calcDistribution();
void tellPoint(int band);
float *distrib = 0;
int *bandCount = 0;
int *distribx100 = 0;
int validDistribution = 0;
Graph distribGraph;
float minDistribScore=1.0, maxDistribScore=1.0;
float minSizeScore=1.0, maxSizeScore=1.0;
float minTolScore=1.0, maxTolScore=1.0;

int Ztrace_clhigh=0;
int tmptest=0;

static int sho;
static void Zgram(), Zsam(), gt55Zap();
void Zrule_menu(), Zclamdump();

static char shufText[14], gapText[14], AText[14], minText[14];
static int shufBox, gapBox, ABox, minBox;
static void scanPz(), pickg9();

char minDistribStr[14], maxDistribStr[14], minSizeStr[14], maxSizeStr[14];
static int txtMin, txtMax;

static void Set_ok() {int i;  
   for (i=0; i<=max_contig; i++) {
       contigs[i].ctgstat=ND; 
   }
   updateproj();
   if (!Zbatch_flag) printf("Complete set all contig status to OK\n");
}
static void Set_nocb() {int i, cnt=0; 
   for (i=0; i<=max_contig; i++) {
      if (contigs[i].seq>0){
       contigs[i].ctgstat=NOCB; 
       cnt++;
      }
   }
   updateproj();
   if (!Zbatch_flag) printf("Complete set %d sequence contig status to NOCB\n", cnt);
}
/****************************************************************
                     DEF: Zrule_menu
***************************************************************/
void Zrule_menu(int who)
{
float col1=1.0, row=1.0;
static MENUOPT amenu[] = {
  { graphDestroy, "Close"},
  { graphPrint,"Print Screen"} ,
  { 0, 0 } };
void CTG_resources();

  sho = who;
  if(graphActivate(grule)){
    graphPop();
  }
  else{
    grule = graphCreate (TEXT_FIT,"The Attic",.2,.2,.35,.35) ;
    graphRegister (PICK, pickg9) ;
  }
  graphClear();
  graphMenu(amenu);

  graphText("Miscellaneous routines:",col1,row);
  row+=1.5;
  if (who==1) {
      graphButton("Unbury selected",Zsel_unbury,2.0,row);
      row+=1.5;
      graphButton("Bury selected",Zsel_bury,2.0,row);
      row+=1.5;
      graphButton("Bury all in selected",Zsel_allbury,2.0,row);
      row+=2.0;
      graphButton("Ctg -> GRAM",Zgram,2.0,row);
      graphButton("Ctg -> SAM",Zsam,15.0,row);
      row+=1.5;
      graphButton("Clone 1 Buried",C1zBur,2.0,row);
      graphButton("Select&Colour",Select_and_Colour,25.0,row);
  }
  else {
      graphButton("Set Sequence Contigs to NoCB",Set_nocb,1.0,row);
      row+=1.5;
      graphButton("Set all Contigs to Okay",Set_ok,1.0,row);
      row+=2.5;
      graphButton("Calculate Distrib",calcDistribution,1.0,row);
      sprintf(minDistribStr, "%f", minDistribScore);
      sprintf(maxDistribStr, "%f", maxDistribScore);
      sprintf(minSizeStr, "%f", minSizeScore);
      sprintf(maxSizeStr, "%f", maxSizeScore);
      txtMin = graphTextEntry(minDistribStr, 4, col1+20, row,scanPz);
      txtMax = graphTextEntry(maxDistribStr, 4, col1+25, row,scanPz);
      txtMin = graphTextEntry(minSizeStr, 4, col1+30, row,scanPz);
      txtMax = graphTextEntry(maxSizeStr, 4, col1+35, row,scanPz);

      row+=1.5;
      graphButton("Histogram to Stdout",histZap,1.0,row);
      row+=1.5;
      graphButton("Greater than 55 bands",gt55Zap,1.0,row);
  }
  sprintf(gapText,"%f", Pz.gap);
  sprintf(AText,"%f", Pz.p_align);
  sprintf(minText,"%d", Pz.minBands);

  row+=2.5;
  graphText("CB parameters:",col1,row);
  row+=2.0; col1++;

  graphText("Percent align",col1,row);
  ABox = graphTextEntry(AText, 5, col1+16, row,scanPz);

  row+=1.5;
  graphText("Maximal gap",col1,row);
  gapBox = graphTextEntry(gapText, 5, col1+16, row,scanPz);

  row+=1.5;
  graphText("Minimal number of bands",col1,row);
  minBox = graphTextEntry(minText, 5, col1+26, row,scanPz);
  row+=1.5;
  graphText("Hit Enter after entering a value",col1,row);

  graphRedraw();
}
/*****************************************************************/
static void scanPz()
{
  if(graphActivate(grule)){
    sscanf(gapText,"%f", &Pz.gap);
    sscanf(AText,"%f", &Pz.p_align);
    sscanf(minText,"%d", &Pz.minBands);
    printf("Gap %f Align %f Min Bands %d\n",Pz.gap, Pz.p_align, Pz.minBands);
  }
}
static void pickg9(int box,double x,double y)
{
  if(box == shufBox) shufBox = graphTextEntry(shufText,8,0,0,NULL);
  else if(box == gapBox) gapBox = graphTextEntry(gapText,8,0,0,NULL);
  else if(box == ABox) ABox = graphTextEntry(AText,14,0,0,NULL);
  else if(box == minBox) minBox = graphTextEntry(minText,14,0,0,NULL);
}
/******************************************************************
		   Get time
******************************************************************/
#include <sys/time.h>
#include <sys/resource.h>
/*#define RUSAGE_SELF 0 */
#define MICRO_UNITS 1000000
struct rusage Start;

void Start_timer()
{
   getrusage(RUSAGE_SELF,&Start);
}

void End_timer()
{
struct rusage info;
long min,sec, min2, sec2;

   getrusage(RUSAGE_SELF,&info);
   min = ((info.ru_utime.tv_sec - Start.ru_utime.tv_sec) / 60);
   sec = ((info.ru_utime.tv_sec  - Start.ru_utime.tv_sec)% 60);

   min2 = ((info.ru_stime.tv_sec - Start.ru_stime.tv_sec) / 60);
   sec2 = ((info.ru_stime.tv_sec  - Start.ru_stime.tv_sec)% 60);

   sprintf(ZBuf,"Time User %ld min %ld sec, System %ld min %ld sec\n",min,sec,
	       min2, sec2);
}
/**************************************************************
		  CTG resources
**************************************************************/
void CTG_resources()
{
struct rusage info;
long min,sec;

   if (getrusage(RUSAGE_SELF,&info)== -1) {
      fprintf(stderr, "Cannot getrusage\n");
      return;
   }
   min = (info.ru_utime.tv_sec / 60);
   sec = (info.ru_utime.tv_sec % 60);
   fprintf(stdout,"Total User  : %ld min, sec %ld, microsec %ld\n",min,sec,
				info.ru_utime.tv_usec);
   min = (info.ru_stime.tv_sec / 60);
   sec = (info.ru_stime.tv_sec % 60);
   fprintf(stdout,"Total System: %ld min, sec %ld, microsec %ld \n",min,sec,
				info.ru_stime.tv_usec);
   fprintf(stdout,"\nPage:  reclaims %ld,  faults %ld,  swaps %ld\n",
	       info.ru_minflt, info.ru_majflt, info.ru_nswap);
   fprintf(stdout,
	 "Resident set size: maximum pages %ld  (= %ld bytes)\n",
	       info.ru_maxrss, info.ru_maxrss * getpagesize());
   fprintf(stdout,"Context switches: voluntary %ld,  involuntary %ld\n",
	     info.ru_nvcsw, info.ru_nivcsw);
   fprintf(stdout,"Shared memory size %ld Swaps %ld\n\n",
	     info.ru_ixrss, info.ru_nswap);
}
/******************************************************************
                   Show time
*******************************************************************/
void Show_time_msg()
{
struct tm *curr_time;
time_t ltime;
char buf[BUFSZ];

     if (time(&ltime) != -1) {
       curr_time = localtime(&ltime);
       if (!strftime(buf, BUFSZ, "%R %a %d %h %Y", curr_time)) {
	    printf("Warning: Not enough room for time stamp.\n");
       }
     }
     else buf[0] = '\0';
     printf("Date: %s\n",buf);
}
/********************************************************************
                   DEF: clamdump - for debugging
***********************************************************************/
void Zclamdump()
{
int i, j;

   printf(
    "\nSet %d N-clone %d N-group %d  n_band %d/%d left %d/%d  right %d/%d\n",
      ZS, Zset[ZS].n_clone, Zset[ZS].n_group, Zset[ZS].n_band, Zset[ZS].max_band,
      Zset[ZS].leftb, Zset[ZS].far_leftb, Zset[ZS].rightb, Zset[ZS].far_rightb);
    printf(" Left: ");
    for (i=0; i<HANG; i++)
        if (Zset[ZS].hang_left[i]!= -1) printf("%d ", Zset[ZS].hang_left[i]);
    printf("\n Right:");
    for (i=0; i<HANG; i++)
        if (Zset[ZS].hang_right[i]!= -1) printf("%d ", Zset[ZS].hang_right[i]);
    printf("\n");

    if (Zset[ZS].n_group == 0) return;

    for (i=Zset[ZS].leftb; i<= Zset[ZS].rightb; i++)
    {
       printf("band %3d avg %4d  high %4d low %4d high-low %3d cnt %2d nocnt %2d clone: ",
         i, Zset[ZS].band[i].avg, Zset[ZS].band[i].high, Zset[ZS].band[i].low,
         Zset[ZS].band[i].high-Zset[ZS].band[i].low, Zset[ZS].band[i].cnt,
         Zset[ZS].band[i].nocnt);
       for (j=0; j < Zset[ZS].band[i].cnt; j++)
            printf("%d ",Zset[ZS].band[i].clone[j]);
       printf("\n");
    }

    for (i=0; i< Zset[ZS].max_band; i++)
       printf("root %3d avg %4d  high %4d low %4d  cnt %4d address %p\n",
        i, Zset[ZS].root[i].avg, Zset[ZS].root[i].high,
         Zset[ZS].root[i].low, Zset[ZS].root[i].cnt, &(Zset[ZS].root[i]));

    for (i=0; i< Zset[ZS].n_group; i++)
       printf("Group %3d left %4d right %4d\n",
        i, Zset[ZS].group[i].leftb, Zset[ZS].group[i].rightb);

    if (Zset[ZS].csort == NULL) return;
    for (i=0; i< Zset[ZS].n_clone; i++) {
         printf("clone %2d: %10s left %2d right %2d extra: ", Zset[ZS].csort[i].i,
           arr(acedata, ZZ.matrix[Zset[ZS].csort[i].i].cin, CLONE).clone,
          ZZ.matrix[Zset[ZS].csort[i].i].leftb,
          ZZ.matrix[Zset[ZS].csort[i].i].rightb);
         for (j=0; j<ZZ.matrix[Zset[ZS].csort[i].i].n_extra; j++)
              printf("%d ", ZZ.matrix[Zset[ZS].csort[i].i].extra[j]);
         printf("\n");
    }
}
/***************************************************************
         fpc->gram
****************************************************************/
static void Zgram()
{
  FILE *fpout;
  CLONE *clone;
  int j;
  char buf[120], buf2[200];
  struct contig *p;

  if(strlen(dirName)>0)
    sprintf(buf, "%s/ctg%d.gram", dirName, currentctg);
  else
    sprintf(buf, "ctg%d.gram", currentctg);

  if((fpout = fopen(buf,"w")) == NULL ){
      printf("file %s does not exist?\n",
             messprintf("%s/%s.fpc", dirName, buf));
      return;
  }
  fprintf (fpout, "GRAM V1.2\n") ;

  for (p=root; p!=NULL; p = p->new)
  {
     clone = arrp(acedata,p->next, CLONE);
     fprintf(fpout, "Clone: %s\n",clone->clone);
     loadCzfast(&C1z, p->next);
     for (j=0; j< C1z.nbands; j++) fprintf(fpout, "%6.4f - - - -\n", (float) C1z.coords[j]/ 1000.0);
  }
  fclose(fpout);
  sprintf(buf2,"Wrote %s",buf);
  displaymess(buf2);
}
/***************************************************************
      fpc->sam
****************************************************************/
static void Zsam()
{
  FILE *fpout;
  CLONE *clone;
  char buf[120], buf2[200];
  struct contig *p;
  struct markertop *topmarkerptr;
  struct marker *markerptr;


  if(strlen(dirName)>0)
    sprintf(buf, "%s/ctg%d.sam", dirName, currentctg);
  else
    sprintf(buf, "ctg%d.sam", currentctg);

  if((fpout = fopen(buf,"w")) == NULL ){
      printf("File %s does not exist?\n",
             messprintf("%s/%s.fpc", dirName, buf));
      return;
  }
  getdatehead(buf2);
  fprintf (fpout, "SAM V2.5 Created: %s\n", buf2) ;

  for (p=root; p!=NULL; p = p->new)
  {
     clone = arrp(acedata,p->next, CLONE);
     if(clone->marker ==NULL) continue;
     fprintf(fpout, "Clone: %s %s\n",clone->clone, clonetype[clone->class]);
     topmarkerptr = clone->marker;
     while(topmarkerptr != NULL){ /*for each marker */
        markerptr = arrp(markerdata,topmarkerptr->markerindex,MARKER);
        if(markerptr != NULL){
          if(topmarkerptr->weak)
            fprintf(fpout,"%16s + L + 1\n", markerptr->marker);
          else
            fprintf(fpout,"%16s + H + 1\n", markerptr->marker);
        }
        topmarkerptr = topmarkerptr->nextmarker;
      }
  }
  fclose(fpout);
  sprintf(buf2,"Wrote %s",buf);
  displaymess(buf2);
}

/**************************************************************
                       histZap
called from Rule Menu from Main analysis
**************************************************************/
void histZap()
{
int v, i, j,  b1, b2;
long int totrange, totunique, cnt;
int tol, totbands, totdbl, totexact;
int min=0;
CLONE *cp;
float avg, adev;
int last, mode, mediam;
int len, cnt_len[20], cnt_bands[20];

struct Zhistogram {
   int cnt, range, tol, peak, dbl, exact;
} *Zhist;

   if (!arrayExists(bands))
        if (fpRead() == -1) return;

   Zhist = (struct Zhistogram *)
             calloc(sizeof(struct Zhistogram), bandmax*2);
   NOMEM2(Zhist, "Zhist");
   totunique=totexact= totrange=totdbl=totbands=0;
   avg=adev=0.0;
   for (i=0; i<bandmax; i++) Zhist[i].exact = Zhist[i].cnt  = Zhist[i].dbl=0;
   for (i=0; i<20; i++)cnt_len[i] = cnt_bands[i]=0;

   for (i=0; i< arrayMax(acedata); i++) {
    cp = arrp(acedata, i, CLONE);
    if (cp->fp==NULL) continue;
    b1 = cp->fp->b1;
    b2 = cp->fp->b2;
    totbands+= b2;
    len = 0;
    last = -1;
    for (j=0; j< b2; j++) {
        v = arr(bands, b1+j-1, int);
        len+=v;
        avg += (float) v;
        Zhist[v].cnt++;
        Zhist[v].tol = tol = Zgtol(v);
        if (abs(v-last) <= tol) {
                Zhist[v].dbl++;
                totdbl++;
                if (v == last) {
                    Zhist[v].exact++;
                    totexact++;
                }
        }
        last = v;
     }
   }
   avg /= totbands;

   for (mode=cnt=0, i=1; i<=bandmax; i++)
   {
       if (Zhist[i].cnt == 0)  continue;

       if (min==0) min=i;
       b1 = MaX(0,i - Zhist[i].tol);
       b2 = MiN(i + Zhist[i].tol, bandmax);
       Zhist[i].range=0;
       for (j=b1; j<=b2; j++) Zhist[i].range += Zhist[j].cnt;
       totunique++;
       totrange+=Zhist[i].range;
       adev += abs(i - avg) * Zhist[i].cnt;
       if (Zhist[i].cnt > Zhist[mode].cnt) mode=i;
   }
   adev /= totbands;

   for (cnt=0, i=1; i<=bandmax && cnt < totbands/2; i++)
        cnt += Zhist[i].cnt;
   mediam = i-1;

   b1=0;
   if (Proj.variable) {
      printf("\nBand  Cnt   Bin   Tol  Dbl   Last\n");
      for (i=1; i<=bandmax; i++)
        if (Zhist[i].cnt > 0) {
             printf("%4d %4d %4d %4d %4d  %3d    %d\n",
                 i, Zhist[i].cnt, Zhist[i].range, Zhist[i].tol,
                 Zhist[i].dbl,
                 i-b1, abs(Zhist[i].cnt - Zhist[i-1].cnt));
             b1=i;
        }
   }
   else {
      b1=0;
      printf("\nBand   Cnt   Bin  Dbl  Last\n");
      for (i=1; i<=bandmax; i++)
        if (Zhist[i].cnt > 0) {
             printf("%4d %4d %4d %4d  %4d %4d\n",
                 i, Zhist[i].cnt, Zhist[i].range,
                   Zhist[i].dbl,i-b1,
                   abs(Zhist[i].cnt - Zhist[i-1].cnt));
             b1=i;
        }
   }
   free(Zhist);
   printf("Total bands %d from %d clones. Average # bands %d.\n",
            totbands, arrayMax(acedata), totbands/arrayMax(acedata));
   printf("Total unique bands %d. Average bin size %d.\n", (int) totunique,
           (int) (totrange/totunique));
   printf("Mode %4d  Minimum %4d  Maximum %4d.\n",mode, min, bandmax);
   printf("Average # doublet %6.3f Total %d\n",
          (float)totdbl/(float)arrayMax(acedata), totdbl);
   printf("Average # exact   %6.3f Total %d\n",
          (float)totexact/(float)arrayMax(acedata), totexact);
   fflush(stdout);
}

void calcDistribution()
{
  int v, i, j,  b1, b2;
  CLONE *cp;
  int numBands = 0;
  int uniqBands = 0;
  unsigned int countDiff[30];
  int maxCount = INT_MIN;

  for (i=0; i<30; i++)
    countDiff[i] = 0;

  if (!arrayExists(bands))
    if (fpRead() == -1)
      return;
  if (distrib == 0) {
    distrib = (float *) calloc(sizeof(float), bandmax*2);
  }
  NOMEM2(distrib, "distrib");

  if (bandCount == 0) {
    bandCount = (int *) calloc(sizeof(int), bandmax*2);
  }
  NOMEM2(bandCount, "bandCount");

  for (i=0; i<=bandmax; i++) {
    distrib[i] = 1;
    bandCount[i] = 0;
  }

  for (i=0; i<arrayMax(acedata); i++) {
    cp = arrp(acedata, i, CLONE);
    if (cp->fp==NULL)
      continue;
    b1 = cp->fp->b1;
    b2 = cp->fp->b2;
    for (j=0; j< b2; j++) {
      v = arr(bands, b1+j-1, int);
      if (v>=0 && v<= bandmax) {
        bandCount[v]++;
        numBands++;
        if (bandCount[v] == 1) {
          uniqBands++;
        }
      }
    }
  }

  for (i=0; i <= bandmax; i++) {
    if (bandCount[i] > maxCount) {
      maxCount = bandCount[i];
    }
    if (bandCount[i] == 0) {
      distrib[i] = 1;
    }
  }

  if (!Zbatch_flag) {
    minDistribScore = atof(minDistribStr);
    maxDistribScore = atof(maxDistribStr);
    minSizeScore = atof(minSizeStr);
    maxSizeScore = atof(maxSizeStr);
  }

  if (Zbatch_flag) {
    fprintf(Zlogfp, "----------------------------------------------\n");
    fprintf(Zlogfp, "Distribution over (%2.1f, %2.1f, %2.1f, %2.1f)\n",
            minDistribScore, maxDistribScore, minSizeScore, maxSizeScore);
  }
  printf("Distribution over (%2.1f, %2.1f, %2.1f, %2.1f)\n",
          minDistribScore, maxDistribScore, minSizeScore, maxSizeScore);

  for (i=0; i <= bandmax; i++) {
    /* distribution score is the minimum plus the range times the percentage
      of the maximum count at that point */
    distrib[i] = minDistribScore + (maxDistribScore-minDistribScore)*
                                   ((float)maxCount - bandCount[i])/maxCount;
    distrib[i] *= minSizeScore + (maxSizeScore-minSizeScore)*
                                   ((float)bandmax-i)/bandmax;
  }

  validDistribution = 1;

  if (!Zbatch_flag) graphArray(&distribGraph, NULL, &bandCount[1], bandmax,
			       "RATE", "COUNT", fileName, "Distribution",
			       tellPoint);
}

void tellPoint(int band)
{
  int v, i, j,  b1, b2;
  CLONE *cp;
  struct list *p=NULL;
  int hitCount=0;

  freelistmem();
  for (i=0; i<arrayMax(acedata); i++) {
    cp = arrp(acedata, i, CLONE);
    if (cp->fp==NULL)
      continue;
    b1 = cp->fp->b1;
    b2 = cp->fp->b2;
    for (j=0; j< b2; j++) {
      v = arr(bands, b1+j-1, int);
      if (v == band) {
        /* this is a clone with the selected band in it. */
        if (listroot == NULL) {
          listroot  = (struct list *)messalloc((sizeof(struct list)));
          p = listroot;
          p->next = NULL;
          p->index = i;
        } else {
          if(p==NULL) p=listroot;
          p->next = (struct list *)messalloc((sizeof(struct list)));
          p=p->next;
          p->next = NULL;
          p->index = i;
        }
        hitCount++;
        continue;
      }
    }
  }

  if (listroot) {
    /* display the results */
    classctg = CLONECLASS;
    displaykeyset(hitCount);
    setDisplayType(5);
    displayProj();
  }
}

static void gt55Zap()
{
CLONE *cp;
int i;

   for (i=0; i< arrayMax(acedata); i++) {
      cp = arrp(acedata, i, CLONE);
      if (cp->fp!=NULL && cp->fp->b2 > 55)
         printf("%s %d\n",cp->clone, cp->fp->b2);
   }
}
